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ABSTRACT 

We present self-consistent triaxial stellar systems that have analytic distribution functions 
(DFs) expressed in terms of the actions. These provide triaxial density profiles with cores or 
cusps at the centre. They are the first self-consistent triaxial models with analytic DFs suitable 
for modelling giant ellipticals and dark haloes. Specifically, we study triaxial models that re¬ 
produce the Hernquist profile from Williams & Evans (2015), as well as flattened isochrones 
of the form proposed by Binney (2014). We explore the kinematics and orbital stmcture 
of these models in some detail. The models typically become more radially anisotropic on 
moving outwards, have velocity ellipsoids aligned in Cartesian coordinates in the centre and 
aligned in spherical polar coordinates in the outer parts. 

In projection, the ellipticity of the isophotes and the position angle of the major axis of 
our models generally changes with radius. So, a natural application is to elliptical galaxies 
that exhibit isophote twisting. As triaxial Stackel models do not show isophote twists, our DFs 
are the first to generate mass density distributions that do exhibit this phenomenon, typically 
with a gradient of « 10°/effective radius, which is comparable to the data. 

Triaxiality is a natural consequence of models that are susceptible to the radial orbit 
instability. We show how a family of spherical models with anisotropy profiles that transition 
from isotropic at the centre to radially anisotropic becomes unstable when the outer anisotropy 
is made sufficiently radial. Models with a larger outer anisotropy can be constructed but are 
found to be triaxial. We argue that the onset of the radial orbit instability can be identified 
with the transition point when adiabatic relaxation yields strongly triaxial rather than weakly 
spherical endpoints. 

Key words: Galaxy, galaxies; kinematics and dynamics - methods; analytical, numerical. 


1 INTRODUCTION 

The case for the importance of triaxiality in galactic dynamics has a 
reasonably long history. Early studies by Binney (1978) and Illing¬ 
worth (1977) proposed that giant elliptical galaxies are generically 
triaxial in shape and this hypothesis was reinforced by the A-body 
simulations of Aarseth & Binney (1978) that showed that aspheri- 
cal initial conditions relaxed to triaxial distributions. More recently, 
triaxiality is believed to be a generic feature of so-called ‘slow ro¬ 
tator’ galaxies (Emsellem et al. 2007; Cappellari et al. 2007; Kra- 
jnovic et al. 2011). 

When tackling the data it becomes difficult to disentangle tri¬ 
axiality from other effects, not least because we are only ever able 
to observe a projection of the distribution. However, even in pro¬ 
jection, triaxiality can give rise to signature effects. In particular, a 
changing axis ratio of concentric ellipsoidal density contours gives 
rise to a twisting of the isophotes when viewed from a general an¬ 
gle. The slow rotators of the SAURON sample (Emsellem et al. 
2007, those with Xu < 0.1) show evidence of isophote twisting. 
In addition to inspecting a single galaxy for evidence of triaxiality, 
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statistics on the shapes of many galaxies yield important informa¬ 
tion. Using the assumption of random viewing angles and a wholly 
axisymmetric population one can recover the distribution of intrin¬ 
sic axis ratios of a sample of galaxies. Eor instance, Tremblay & 
MeiTitt (1995) showed from a sample of 171 bright ellipticals that 
wholly oblate or prolate populations were ruled out and a triaxial 
population favoured. More recently, Weijmans et al. (2014) per¬ 
formed a similar procedure on the ATLAS-3D sample (Cappellari 
et al. 2011) and showed there were strong indications of triaxiality 
in the sample of slow rotators. 

The recent SAURON and ATLAS-3D projects have com¬ 
plemented photometry with line-of-sight velocity distributions for 
samples of nearby galaxies. This additional information helps to 
unravel the internal dynamics, and hence the intrinsic shapes, of 
the galaxies. Binney (1985) proposed that rotation along the minor 
axis of projected elliptical density contours or instead kinematic 
misalignment (misalignment of the minor axis and the angular mo¬ 
mentum vector) is indicative of triaxiality. Eranx et al. (1991) in¬ 
spected the statistics of minor-axis rotation for a sample of galax¬ 
ies and found 35 per cent of their sample had kinematic mis¬ 
alignment. Weijmans et al. (2014) used kinematic misalignment to 
separate off potential triaxial candidates in their sample. 
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So far our discussion of the evidence for triaxiality has been 
limited to the visible component of galaxies. The shapes of dark 
matter haloes remains an interesting open question. Many A'^-body 
simulators (Jing & Suto 2002; Bailin & Steinmetz 2005; Allgood 
et al. 2006) have demonstrated the range of possible triaxial dark 
matter haloes in cosmological simulations. However, the models 
relax to axisymmetry or sphericity when baryons are included in 
the simulations (Udry 1993; Dubinski 1994; Barnes & Hernquist 
1996; Kazantzidis et al. 2004; Debattista et al. 2008) or when a 
central black hole forms (Merritt & Quinlan 1998). The case for 
triaxiality of the Milky Way’s dark matter halo is motivated by stud¬ 
ies of the Sagittarius stream (Law & Majewski 2010; Vera-Ciro & 
Helmi 2013; Deg & Widrow 2013). However, several of these re¬ 
sults are unsettling in that the short-axis of the dark matter halo lies 
in the disc plane (Debattista et al. 2013). Unhappily, the Sagittarius 
stream has proved hard to model, with no single fit able to explain 
the wealth of new data. 

1.1 Dynamical models 

Although there is considerable evidence for the presence of triaxi¬ 
ality in galaxies, the range of modelling tools available is severely 
limited. The two greatest developments in the understanding of tri¬ 
axial stellar systems were the work of de Zeeuw (1985) on Stackel 
or separable potentials and Schwarzschild (1979) on the construc¬ 
tion of numerical dynamical equilibria, de Zeeuw (1985) showed 
that the only Stackel potential in which the isodensity surfaces are 
concentric ellipsoids is the perfect ellipsoid and provided key in¬ 
sights into the classes of orbits that are expected to arise in generic 
triaxial potentials. Schwarzschild (1979) showed how a given tri¬ 
axial density-potential pair could be self-consistently generated 
by linear superposition of numerically integrated orbits. Statler 
(1987) used this method to construct numerical dynamical mod¬ 
els of de Zeeuw’s perfect ellipsoid and demonstrated that they 
genetically produce kinematic misalignment when viewed in pro¬ 
jection. However, as shown by Franx (1988), triaxial Stackel mod¬ 
els are unusual in that they are unable to produce isophote twisting. 
Schwarzschild’s method has proved useful in modelling external 
galaxies, but the only triaxial fit to an elliptical galaxy remains the 
work of van den Bosch et al. (2008) on NGC 4365 - though the 
related problem of modelling rotating triaxial bars has seen some 
progress (Zhao 1996; Hafner et al. 2000; Wang et al. 2013). Triaxial 
models have also be constructed using a made-to-measure approach 
(Syer & Tremaine 1996; Dehnen 2009), though here the only fits 
to data are for the Milky Way bar (Long et al. 2013; Portail et al. 
2015). 

An alternative approach to the construction of dynamical equi¬ 
libria is through the use of integrals of motion. The Jeans theo¬ 
rem states that an equilibrium distribution function (DF) is solely a 
function of the integrals of motion. In spherical systems the classi¬ 
cal integrals of the energy, E, and angular momentum, L, can be 
used to construct equilibrium models. For instance, Eddington in¬ 
version allows us to construct the isotropic f{E) model given any 
radial density profile p(r) and these calculations may incorporate 
anisotropy in a variety of ways (Osipkov 1979; Merritt 1985; Evans 
& An 2006; Binney & Tremaine 2008). Moving from sphericity to 
axisymmetry, it has long been realised that, in addition to the en¬ 
ergy and a single-component of the angular momentum (e.g. the 
2 -component Lz), numerically integrated orbits obey a third non- 
classical integral, I 3 (e.g., Ollongren 1962) and such a dependence 
in the DF seems necessary to reproduce the dynamics of the Milky 
Way and external galaxies (e.g., Binney 1976; Davies et al. 1983). 


In general potentials, a global analytic third integral does not exist, 
but in axisymmetric Stackel potentials this integral can be written 
down (e.g., Lynden-Bell 1962; de Zeeuw 1985). Similarly, in triax¬ 
ial Stackel potentials, we can find a second non-classical integral I 2 
(a generalisation of a component of the angular momentum) to give 
three globally defined analytic integrals of motion. However, no an¬ 
alytic DFs based on the three integrals have ever been constructed 
for triaxial Stackel models, with the exception of the very idealized 
case when the tube orbits are all infinitesimally thin (Hunter & de 
Zeeuw 1992). There are some DFs known for rotating triaxial stel¬ 
lar systems, though only for rather unrealistic density distributions 
such as homogeneous ellipsoids or polytropes (Ereeman 1966; Van- 
dervoort 1980) 

1.2 Action-based distribution functions 

As any function of the integrals of motion is also an integral of mo¬ 
tion, we are free to take any choice of these functions as the basis 
for our DF. A natural choice is the action coordinates, J. Along 
with the angle variables, they form a canonical set of coordinates. 
They possess the following advantages over the classical integrals: 
the range of each action is independent of the other actions, they are 
adiabatic invariants (crucial for the application in this paper) and 
they possess a physical meaning e.g. the ‘radial’ action Jr describes 
the extent of radial excursions (see Binney & Tremaine 2008, for 
more information). As the actions are adiabatic invariants, we are 
able to propose a form for the DF and iteratively converge to a self- 
consistent solution. However, in order to do this we must evalu¬ 
ate the properties of a given f{J) model so we require routines to 
find J given {x,v). In spherical potentials the actions are given 
by J = {Jr, J(fi, Je) = {Jr, Lz,L — \Lz\)- Jr can be calculated 
as a quadrature for any spherical potential. In recent years, several 
algorithms have been developed to estimate the actions in general 
potentials. Sanders (2012) proposed an algorithm where an axisym¬ 
metric potential was locally fitted by a Stackel potential and then 
the actions estimated as those in the best-fitting Stackel potential. 
A similar routine was presented by Binney (2012) who rewrote the 
equations for the actions in a Stackel potential in such a way that 
they could be applied to a general potential. The accuracy of both 
of these approaches relied on the target potential being sufficiently 
close to a Stackel potential but crucially only locally. Sanders & 
Binney (2015) built on the work of Binney (2012) by generalis¬ 
ing the method to triaxial potentials. In addition to these approx¬ 
imate methods for finding the actions, Sanders & Binney (2014) 
presented a more accurate approach for general triaxial potentials 
that relied on the construction of a generating function from time 
samples of an orbit integration. 

Given the tools now available for the calculation of actions, 
it is natural to begin constructing self-consistent action-based DFs. 
Binney (2014) demonstrated how the isochrone DF could be flat¬ 
tened to create a family of axisymmetric action-based models. 
More recently, both Posti et al. (2015) and Williams & Evans 
(2015) have proposed families of action-based distribution func¬ 
tions that reproduced the density profiles and kinematics of a range 
of spherical potentials. With the tools in hand to estimate actions in 
triaxial potentials, this paper makes the obvious next step of gener¬ 
alizing these models to triaxiality with the aim of producing models 
that are appropriate for the modelling of the triaxial structures dis¬ 
cussed in the previous sections. We demonstrate how triaxial self- 
consistent models can be generated and inspect the properties of a 
few illustrative cases. 

In Section 2, we detail all the pieces of machinery required 


(gl 2015 RAS, MNRAS 000, 1-16 


Self-consistent triaxial models 3 


to construct self-consistent triaxial models. Some technical details 
are reserved for three appendices. In Section 3, we present triax¬ 
ial versions of two recent action-based distribution functions in the 
literature. In Section 4, we discuss two novelties of the presented 
models: isophote twisting and the ability to explore the radial sta¬ 
bility of the models. Finally, we discuss the importance of reso¬ 
nances and chaos for these models and other uses of the models in 
Section 5 and draw our conclusions in Section 6 . 


2 CONSTRUCTING A SELF-CONSISTENT TRIAXIAL 
MODEL 

In order to find the self-consistent potential 4? for a given f{J) we 
require a considerable amount of machinery. Here we will outline 
and review the various pieces of machinery required. 


2.1 Review of the action estimation schemes 

For the calculation of properties of a given f{J), we require algo¬ 
rithms for the computation of J from (x, v) given some general 
potential. Actions only exist for integrable potentials and chaos can 
play a significant role in many triaxial potentials (Merritt & Frid¬ 
man 1996; Siopis & Kandrup 2000). We will assume in what fol¬ 
lows that all orbits can be labelled by an action (or some approxi¬ 
mate action) and defer discussion of the presence of resonances and 
chaos to section 5. 

We will use two algorithms for the computation of actions. 
The first is the triaxial Stackel fudge method introduced in Sanders 
& Binney (2015) and is appropriate for the rapid generation of 
models. The second is a slower but more accurate method using 
the generating function from a analytic set of angle-actions to the 
target calculated from orbit integration. This second method was 
detailed in Sanders & Binney (2014) and will be useful for check¬ 
ing the results of our calculations using the Stackel fudge. Here, 
we give an outline of the first of these methods and detail any al¬ 
terations made for the current application. The second method is 
detailed in Appendix C where a cross-check of our calculations is 
presented. 


2.2 Triaxial Stackel fudge 

Sanders & Binney (2015) introduced an algorithm for rapidly esti¬ 
mating the actions in a general triaxial potential. It built on the work 
of Binney (2012) and operates by assuming the target potential is 
sufficiently close to a Stackel potential over the region a given orbit 
explores. 

The Hamilton-Jacobi equations are separable in a Stackel po¬ 
tential such that the actions can be expressed as ID quadratures. 
The canonical coordinates in which the equations separate are 
{t,Pt) where r are ellipsoidal coordinates (A,/r, given by the 
roots of 


T + a 


-F 


y 


= 1 . 


( 1 ) 


T + p r -F 7 

Here, {x, y, z) are Cartesian coordinates and a, /3 ,7 are parameters 
that describe the confocal coordinate system. We impose a < /? < 
7 such that X is the long axis and 2 : the short axis. 

The most general Stackel potential, "Fs, can be expressed as 




/(A) 


(X/iv) 


(A-/r)(t/-A)’ 


( 2 ) 


where {Xyu) denotes cyclic permutations of the three variables. 
<I >5 is composed of a single function of one variable, /(r). The 
equations of motion can be written as 

2(r -F a)(T -F /3)(t -F 7)Pr = — ra -F fe -F /(r) (3) 

where E is the energy and a and b are separation constants, pr is 
solely a function of r such that the actions are given by ID quadra¬ 
tures 

Jt = - [ dr |pr('r)|. (4) 

7t 


(r_, r+) are the roots of the right hand side of equation (3). Note 
that as in Sanders & Binney (2015) we have defined a full oscilla¬ 
tion in r to be two full oscillations from r_ to r+ and back again. 
This gives twice the true radial action for loop orbits but means 
orbits continuously fill action space (Binney & Spergel 1984). 
Given a general triaxial potential, we define the quantities 

Xx{X,y,v) = - p){v - p,u), 

Xm(A, >x) = (p- u){X - p)^{X, p, v), (5) 

Xi,{X,p,v) = {v - X){p - v)^{X, p,v). 


If 4? were a Stackel potential, these quantities would be given by, 
for instance. 


Xx{X,p,u) 


/(A) - 


Therefore, for a general potential, we can write 


( 6 ) 


/(t) ~ XTiX,p,l')-\-C-rr-\-Dr, 


(7) 


where Cr and Dr are constants provided we always evaluate Xt 
with two of the ellipsoidal coordinates fixed. For instance, we al¬ 
ways evaluate xa at fixed p and u. 

When we substitute these expressions into equation (3), we 

find 


2(r-Fa)(r-F/3)(r-F7)pT = E - xAr + Br EXriX, p,v). (8) 


For each r coordinate, there are two new integrals of motion given 
by Ar = a — Cr and Br = b + Dr - Using a single 6 D coordinate 
and a choice of coordinate system gives us a single constraint on 
a combination of Ar Br ■ Due to the separability of the Hamilton- 
Jacobi equation, the partial derivative of the Hamiltonian with re¬ 
spect to any of the ellipsoidal coordinates is zero for a true Stackel 
potential. Setting it equal to zero for a general potential gives a 
further constraint on At and Br allowing us to solve for these in¬ 
tegrals given only a single {x,v) coordinate. We have therefore 
produced approximate pr (r) equations of motion which may be 
integrated to estimate the actions. 

In the above algorithm, the only parameters we can control 
are a and /3 which determine the location of the foci. We set 7 = 
— 1 without any consequences. These are chosen on an orbit-by- 
orbit basis so we are free to change them every time we require 
another action. We give details of the method employed to do this 
in Appendix A. 


2.3 Adiabatic relaxation 

A DF must obey the collisionless Boltzmann equation. For some 
applications such as the stellar halo, we may consider the DF of 
a tracer population with negligible mass that lives in the potential 
generated by a more dominant component. However, when treating 
massive components of galaxies, it is preferable for the model to 
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self-consistently generate its own potential. This poses a problem 
for action-based DFs. In order to find the density and the poten¬ 
tial of an action-based DF, we must know the potential in order to 
calculate the actions. Therefore, the models must be constructed 
iteratively. 

As the actions are adiabatic invariants, the form of the DF does 
not change under slow changes of the potential. We therefore pose 
an initial guess of the potential "To and find the potential $1 of 
the DF in $ 0 . We repeat this procedure with T'o replaced by <l?i 
until the difference between the two potentials is less than ~ one 
per cent everywhere. This procedure is not guaranteed to converge 
if >l>o is very different from <l?i. To aid with convergence Binney 
(2014) used a linear combination of $0 and $1 as the next trial po¬ 
tential. We did not adopt that procedure here but found our models 
converged within ~ 8 iterations. 

2.4 Multipole expansion 

At each stage in the iterative procedure, we require the potential 
calculated from the density. Here, we are working with general tri- 
axial mass distributions so we evaluate the potential using a multi¬ 
pole expansion (Binney & Tremaine 2008). We discuss the details 
of this in Appendix B. 

2.5 Tests 

The code was tested by reproducing spherical versions of the mod¬ 
els detailed in the next section allowing for triaxiality and compar¬ 
ing to the models constructed by limiting the models to sphericity 
(i.e. using only I — 0 and m = 0 terms in the multipole expan¬ 
sion and calculating the actions in a spherical potential). Tests of 
the multipole expansion code are given in Appendix B. The action¬ 
finding code has been tested in Sanders & Binney (2015) and it 
was found that the actions were accurate to < 10 per cent for a typ¬ 
ical box orbit and < 5 per cent for a typical loop orbit in a triaxial 
NFW potential. A comparison of the density calculation with two 
different action estimation approaches is presented in Appendix C. 


3 ACTION-BASED MODELS 

Due to the limited number of cases for which analytic self- 
consistent action-based DFs can be computed the study of action- 
based DFs has been fairly limited. In this section, we give details of 
several recent models which will be of use in constructing triaxial 
models. We consider DFs that are functions of the absolute values 
of the actions. Therefore, these models do not have a streaming 
velocity. This is not a necessary requirement for the models but it 
simplifies their construction. In all formulae, M is the mass of the 
model and A/” is the normalization given by 

JV = M-\2nf f d®J/(J) 

(9) 

COO coo coo ' 

= M-\2nfT dJr dJ^ dJgf{J). 

Jo Jo Jo 

The definition of the actions in the triaxial case (i.e. the radial action 
for the loop orbits is multiplied by two) means the integral includes 
loop orbits rotating in both senses so = 1 but in the spherical 
case we must multiply the integral by a factor of = 2 to include 
orbits rotating both clockwise and anticlockwise. We introduce a 
change of variables Si = Ji/{Jo + Ji) to map the infinite limit 
to (0,1) and calculate the integral numerically using the Divonne 


routine in the CUBA package (Hahn 2005) (all integrals over the DF 
are performed with this package). Additionally in the triaxial case 
we evaluate the distribution function as f(^Jr,J<f,,Jg) such that 
we correctly map onto the spherical case. 

Sanders & Binney (2015) inspected a simple DF given by 

fsB{J) = MJV{Jo+JZ{J)r, (10) 

where 

JZ{J) = Jr CLip\J,p\ O-elJo]- (11) 

Here, the coefficients ai are constants and Jo is a scale action. The 
models were evaluated in a fixed Navarro-Frenk-White potential 
and so were not self-consistent. They have a density core and den¬ 
sity fall-off like r~^ and the coefficients controlled the kinematics. 
Two models were studied and it was demonstrated that the Jeans 
equation was satisfied to good accuracy. These models are simple, 
but lack flexibility. Here we examine two other models from the 
literature which we will study further with our new machinery. 


3.1 Williams & Evans models 


Williams & Evans (2015) proposed a family of action-based mod¬ 
els with double power-law density profiles of the form 

p oc r“(ro J-(12) 


and tunable anisotropy profiles. The models built on the work of 
Williams et al. (2014) who demonstrated that in scale-free spherical 
potentials (3> oc r"), the Hamiltonian is well approximated by a 
homogeneous function of degree one in the actions i.e. H{J) oc 
(L -I- Dweb Jr)^ where ( — 2e/{e -I- 2). Dvveb is a coefficient 
that can be calculated by considering the limiting cases Jr = 0 and 
L — 0. Adjusting Dweb away from the isotropic value produces 
models with constant non-zero anisotropy. 

For a double power law in the density, the potential in the near 
and far field is well approximated by a scale-free power law such 
that an appropriate isotropic DF can be constructed by stitching to¬ 
gether two scale-free H{J) as 


/we(J) =AfM 


Ti\J\)£{jy 


ji + £{jy] 


(ri-5)/2 ’ 


(13) 


where 


JZ{J) = L + D{\J\)Jr.. (14) 

Here, Jo is a scale action related approximately to the scale radius 
To of the model as Jo = y/GMro and 5 = (6 — a)/(4 — a) and 
p = 26 — 3. Both D(| J|) and r(| J|) are transit functions of the 
form 

^<'•'11 = <‘S 

where | J| = v/J^”+TJ. Td J| ) evolves from To to 7i over a scale 
Jr- r(| J|) controls the weight of the near-field part of the model 
relative to the far-field part. It acts to adjust the location of the 
break-radius in the density profile. D(| J|) controls the anisotropy 
of the model. Note that the scale-action is summed in quadrature 
with C. This was found by Williams & Evans (2015) to produce 
better fits to the curvature of the density profile around the break 
radius. The model presented by Williams & Evans (2015) is very 
similar to that presented by Posti et al. (2015). However, the model 
of Williams & Evans (2015) is more flexible in the range of density 
profiles and kinematics that can be fitted so we adopt it here. 
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In the spherical case the actions and Jg are combined into 
L = I + |>/e| so Jg = L — I J<^|. To generalize these models 
beyond sphericity, we set 

D{\J\)\Jr\ + \Jf + \Je\/q.. (16) 

qz acts to flatten the model in the 2 -direction. Some combinations 
of parameters will produce axisymmetric models, whilst others will 
be triaxial, but there is no way of finding the symmetry of the model 
without a full calculation. 


3.1.1 Example model 


As a demonstration of the methodology, we will construct triaxial 
analogues of the Hemquist (1990) profile using equation (13). We 
refer to these as WEH models, as they use the ansatz introduced by 
Williams & Evans (2015). 

A spherical Hemquist profile has the density-potential pair 


PH(r) 

$H(r) 


M 1 
27t r(ro -f rf ’ 
_ GM 
r -I- ro' 


(17) 


We set G = M = ro = 1. To constmct an appropriate action- 
based model for this profile, we take the parameters from Williams 
& Evans (2015) as a starting point. They find an isotropic Hemquist 
profile is given by <5 = 5/3, rj — 5, Do = Dweb = 1.814, 
Hi = 1, To = 0.378, Ti = 1, Jo = Jt = 1 and Jd = 0.41. 
We then introduce a flattening qz = 0.4. This simple alteration 
produces a model that is oblate in the far-field but weakly prolate 
in the centre. As our action-finding algorithm is not designed to 
find the actions in prolate cases, these results are perhaps not to 
be completely tmsted. To encourage the model to be oblate at the 
centre, we increased the central radial anisotropy by decreasing Dq. 
As noted in Williams & Evans (2015) when altering Di, Ti must 
be adjusted accordingly to retain the required density profile as 


To ^To( 
Ti Ti 


1 + Hwbb N 
1 + Ho / 
2 \-v 

1 + dJ ■ 


-s 


( 18 ) 


For Dq = 1, the model is triaxial at the centre (we will dis¬ 
cuss under what conditions the models become triaxial later). We 
use as our initial potential guess a Hemquist profile flattened in the 
potential by b/a = 0.98 and c/a = 0.95. In Fig. I, we show the 
convergence of the Dq = 1 model towards self-consistency. Note 
also that at each radial point, there is a spread of points correspond¬ 
ing to a range of spherical polar angles. 

In Figure 2, we plot the density contours of the model in the 
{x, y) plane and (x, z) plane along with the axis ratios of the best¬ 
fitting ellipses to the density and potential contours, b/a is the ra¬ 
tio in the (x, y) plane, whilst c/a is the ratio in the (x, z) plane 
(note here we could simply use the ratio of the intercepts of the 
coordinate axes, but as the models exhibit a slight peanut shape 
this gives a false impression of the shape of the contours). We 
can see that the model is triaxial at the centre with central values 
( 6 /o)p(r = ro/10) « 0.7 and {c/a)p{r = ro/10) « 0.4. For 
r > lOro, the potential is much more spherical, whilst the den¬ 
sity is still flattened. At these large radii, there is little mass. As 
the total mass is finite, the behaviour at large radii mimics that of 
a tracer population in a monopole potential. The density contours 
in (x, z) are not elliptically shaped, but take on a slight peanut 
shape. Emsellem et al. (2011) shows that the slow rotators in the 



10 -^ 10 -‘ 10 “ 10 * 10 ^ 
r 


Figure 1. Convergence and final profiles of a triaxial WEH model: the top 
panel shows three iterations of the density calculation at arbitrary multi¬ 
plicative offsets for visibility (the final iteration has no offset) and the bot¬ 
tom panel shows the corresponding potentials with arbitrary additive off¬ 
sets. The blue dots show the calculation in the initial Hemquist potential 
(shown in grey), the green down-facing carets show the first iteration, red 
crosses the third and final iteration. Note we show the density calculated at 
several {cj>, 6) values for each spherical radius r to give an indication of the 
triaxiality of the model. The positions are in units of ro, the density in units 
of M/xq and the potential in units of GM/vq. 


ATLAS-3D sample (i.e. possible triaxial galaxies) have a boxiness 
ratio ai/a consistent with pure elliptical or slightly boxy. Those 
that have slightly discy contours appear to have a kinematically de¬ 
coupled core. It is therefore interesting that peanut-shaped distribu¬ 
tions arise in the models (and also in the models of Schwarzschild 
(1979)), but not in nature. However, we should note that the de¬ 
gree of triaxiality presented here is larger than many galaxies are 
believed to have. Also, we will see later that in projection these 
models look elliptical. 

In Figure 3, we show the density profile along a radial ray at 
the spherical polar angles </ = f and 0 = j decomposed in terms 
of the orbit classification. We see that in the centre of the model, the 
majority of the density is contributed by the box orbits whilst in the 
outskirts the long and short axis loops contribute equally. The or¬ 
bits are classified based on their limits in the Stackel fudge scheme 
as described in Sanders & Binney (2015). A small fraction of or¬ 
bits have boundaries that do not correspond to any orbital class in a 
Stackel potential so we label them as approximately box, short-axis 
loop or long-axis loop depending on the relative magnitude of their 
actions. Our orbit classification method artificially imposes regu¬ 
larity on the orbital contribution to the model as some orbits will 
be chaotic but incorrectly given a regular label. There exist several 
methods for ascertaining whether an orbit is regular or chaotic (e.g. 
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Si 




Figure 2. Density and shape of the triaxial WEH model: the left panel shows a slice of the density in the {x, y) plane and the middle panel shows the equivalent 
in the (a;, z) plane. The blue colours indicate the base 10 logarithm of the density relative to the central value. The insets show zoom-ins for |a;| < 1, \y\ < 1 
and |a;| < 1,|2| < 1. In the right-hand panel, we show the axis ratios of the fitted ellipses in these two planes with solid lines [6/a corresponding to (x, y) 
and c/a (x, z)]. The dashed lines show the equivalent for the potential. The positions are in units of ro and the density is in units of M/xq. 



Figure 3. Contributions to the density in the triaxial WEH model along a 
radial ray with (f> = 0 = The boxes contribute most at the centre, 

whilst the loop orbits contribute equally in the outskirts where the potential 
is roughly spherical. Also shown is a true Hernquist profile offset by a factor 
of three for visibility. The positions are in units of ro and the density is in 
units of M/ro. 


Lyapunov exponents - see Vasiliev (2013a) for a nice summary or 
the small alignment index [SALI] Skokos et al. (2003) that was ap¬ 
plied recently to characterise the orhits in an axisymmetric galactic 
potential by Zotos (2014)) but these require long orbit integration. 
It seems for non-rotating galactic potentials of interest the majority 
of weight is on orbits that appear regular over a Hubble time (see 
Section 5 for more discussion) such that our orbital decomposition 
is meaningful. 

In Figure 4, we show the velocity dispersions in the (x, y) 
and (x, z) planes, and ayy produce approximately elliptical 
contours whilst produces a narrow waisted distribution in the 
(*, z) plane with the dispersion falling off much more slowly along 
the 2 axis. We also plot the velocity ellipses in the same two planes. 
Clearly the velocity ellipses are very radial for 2 > 0. In the (x, z) 


plane, the ellipses are near radially-aligned for most of the explored 
space. At large radii, the ellipses become more radially aligned and 
at small radii they appear to become aligned with the Cartesian co¬ 
ordinates. In the (a;, y) plane, we have a similar situation but the 
ellipses are near circular so any alignment of the ellipses in this 
plane is subtle. This configuration is seen in models with Stackel 
potentials as the velocity ellipsoids are aligned with the confocal 
ellipsoidal coordinate system. Fig. 3 shows that at the centre of the 
model the box orbits dominate so the coordinate system is naturally 
aligned with Cartesian coordinates whilst beyond the scale radius 
loop orbits become more significant and the velocity ellipses nat¬ 
urally become more spherical aligned. Although we have used a 
Stackel-based approach to estimate the actions, we have allowed 
the focal distances of the coordinate system to be a function of en¬ 
ergy so we believe that this configuration is not a consequence of 
our calculation but a genuine feature of the models 


3.2 Binney’s flattened isochrones 

The spherical isochrone density-potential pair is given by 


Pi(r) 

$i(r) 


M 3(ro -b ra)rl - r^(ro -|- Srg) 
47t ri(ro + ra)3 

GM 

ro+Xa’ 


(19) 


where Va = ^/xq + r'^ (Binney & Tremaine 2008). The isochrone 
density profile has a constant core and falls off as at large radii. 
The isochrone potential is convenient as it is the most general po¬ 
tential in which the actions can be analytically calculated (Evans 
et al. 1990). The Hamiltonian as a function of the actions is 


H'i(J) = 


{GMf 


( 20 ) 
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Figure 4. Velocity dispersions and ellipses of triaxial WEH model: in the top row the first panel shows in the {x, y) plane, the second panel shows <jyy 
in the (x, y) plane, the third panel shows in the {x, z) plane and the fourth panel shows in the (a:, z) plane. In the second row, the first panel shows 
the velocity ellipses in the (x, y) plane, the second panel shows the tilt with respect to the radial direction (drawn in light grey in the first panel) in the (x, y) 
plane, the third panel shows the velocity ellipses in the (x, z) plane and the fourth panel shows shows the tilt with respect to the radial in the (a:, z) plane. 
The colours show the magnitude of the velocity dispersion along the major axis of each ellipse and only the shape of the ellipse is important. The sizes of the 
ellipses are arbitrary. The bottom row shows a zoom-in of the second row. The positions are in units of ro and the dispersions are in units of sjGM/ vq. 


Henon (1960) demonstrated via Eddington inversion that the 
isotropic distribution function is given by 


hm 


_V_ Vn 

^(27t)3(GMro)3/2 [2(1 - H)]‘^ 


[27 - 66 U + 320U'^ 


constants a,j>o and ago - The coefficients are given by 
ar{J) = (1 - ip)ao + iparo, 
a^( J) = (1 - ip)ao + 
ag(J) = ago, 


(24) 


- 240-H^ + 647f'‘ + 3(16-H^ + 28-H-9)-^ 


- -H) 


]■ 


where 


H = - 


Hiro 
GM ■ 


where 


(21) 

ao{J) 


(«^) 

(22) 

V'(J) 


fir 


{ago — 1), 


(11 + fl 

— -^{a<l >0 + Olgo 

v>( J) = tanh(J/VGMro). 


(25) 


Note here A/” = 1 but we retain it as it is necessary when adjusting 
the DF. Binney (2014) proposed a scheme to flatten these models 
such that axisymmetric analogues of the isochrone could be con¬ 
structed. The density of the model in planes of constant energy is 
adjusted by scaling each of the actions by The new model is 
given by 


In these expressions, the frequencies, Q-i, are evaluated in the spher¬ 
ical isochrone potential at the actions J = (J, J, J) which is the 
barycentre of a plane of constant energy given by^ 


^ 2GM 



(26) 


jYi{J^ — ara^pag J\{aj-JT, agjg{. (23) 

In order to retain the spherically-averaged density profile of the 
isochrone, the Qi are not independent. We are free to choose two 


3 Note there is a typo in equation (11) of Binney (2014) such that the 
expression inside the squai'e root should read ^ + SGMro, not 

- 3GMro. 
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In this equation E is evaluated using H{J) for the spherical 
isochrone. Expressions for the frequencies are given in Binney & 
Tremaine (2008). The linear combination in equation (24) was in¬ 
troduced by Binney (2014) to produce models that did not become 
prolate at the centre. In this way ar tends to at low barycentric 
action. 

3.2.1 Example model 

We now construct a triaxial generalization of one of these models. 
In order to induce triaxiality, we require ago > 1 such that the 
model is flattened in the 2 : direction and aro < o <#>0 such that the 
model is radially biased. However, we also require a,;j,o + «eo < 3 
such that Qro > 0 everywhere. 

We initially tried the a,^o = l,ago = 1.4 model presented 
in Binney (2014) but this produced a near-axisymmetric model. 
We instead opted for ago = 1.3 and a,^o = 1-6 and we set 
G — M = ro = 1. This produced a weakly triaxial model. In 
Figure 5, we plot the density contours of the model in the (x, y) 
plane and {x, z) plane along with the axis ratios of the best-fitting 
ellipses to the density and potential contours. We can see that 
model is triaxial at the centre with {b/a)p(r = ro/ 10 ) ~ 0.8 and 
{c/a)p{r = ro/10) « 0.7. For r > lOro the potential is tending 
towards spherical whilst the density is still flattened. We see the 
strong core in the density contours of the model and the slight boxy 
shape of the density contours in the (x, z) plane. 

In Figure 6 , we show the density profile along a radial line 
at the spherical polar angles <j) — j and 0 — ^ decomposed in 
terms of the orbit classification. At the centre, nearly all the den¬ 
sity is contributed by the box orbits. The contribution to the density 
from the loop orbits turns over around r = ro due the harmonic 
core of the potential. For ro < r < 3ro the short-axis loop orbits 
and the box orbits are contributing equally whilst beyond r = 3ro 
the long-axis loop orbits dominate. In the outskirts, the long-axis 
loops contribute most with very weak contributions from the box 
orbits and the short-axis loops. This combination is clearly neces¬ 
sary to produce the observed triaxial shape in the outskirts. It is 
interesting to note the similarities and differences from the Hem- 
quist model. In both cases the density of box orbits outweights the 
loop orbits by two orders of magnitude. However, for the isochrone 
model the long and short axis loops contribute equally whilst for the 
Hemquist model the long axis loop orbit contributes over an order 
of magnitude less density than the short axis loop orbits. At large 
radii the box orbits contribute weakly in both cases with the loop 
orbits contributing most to the density. However, in the Hemquist 
model the short and long axis loop orbits are contributing equally at 
large radii whilst in the isochrone model the long-axis loop orbits 
are dominating the density. Note that in this case, unlike with the 
WEH model, the total density profile differs from an isochrone pro¬ 
file. The models will only exactly match for small |ai — 1| which is 
not satisfied for this model as q_r « 3 — — Qz « 0.1. However, 

the inner and outer slopes match well. 

In Figure 7, we show the velocity dispersions in the (x, y) 
and (x, z) planes, and ayy produce approximately elliptical 
contours within r « 2ro as with the Hemquist case. However, the 
velocity dispersions at the centre do not rise nearly as steeply as the 
in the Hemquist case. Interestingly, the contours break from el- 
lipticity for y > 3ro as the outermost contour shown has a positive 
curvature. Also ayy seems to have a narrow waisted distribution for 
X > 3ro. As in the Hemquist case, a\^ produces a narrow-waisted 
distribution in the {x, z) plane with the dispersion falling off much 
more slowly along the 2 axis. We also plot the velocity ellipses in 



Figure 6. Contributions to the density in the triaxial Binney (2014) 
isochrone model along a radial line with <j> = 8 = Also shown 

is an isochrone profile offset by a factor of three for visibility. The positions 
are in units of ro and the density is in units of M/ro. 


the same two planes. Clearly the velocity ellipses are very radial in 
both planes. The stmcture of the radial alignment in both planes is 
similar to the Hemquist case with the largest deviation from radial 
alignment occurring at a; = 0 and small y and 2 . 


4 APPLICATIONS OF MODELS 

We present two applications of our new triaxial models. 


4.1 Isophote twisting 

One of the main lines of evidence for the triaxiality of some el¬ 
liptical galaxies is the observation of isophote twisting (Bertola 
1981; Emsellem et al. 2007). Isophote twisting is the change in 
the position angle of the major axis of the isophotal contours. In 
fast rotating galaxies or bulges (e.g.. Stark 1977), this could be 
attributed to an inner bar/disc structure, whereas in slow rotators, 
the phenomenon is often attributed to triaxiality (e.g. Williams & 
Schwarzschild 1979). The projection of a triaxial model with el¬ 
lipsoidal axis ratios varying with distance from the centre natu¬ 
rally produces gradual isophotal twists. Another explanation for the 
isophote twisting is an intrinsic twisting possibly due to a recent 
interaction. Several authors have studied the statistics of isophote 
twisting of large samples of galaxies (Carter 1978; Benacchio & 
Galletta 1980; Leach 1981). Fasano & Bonoli (1989) analysed a 
sample of galaxies and concluded that tidal interactions only played 
a small role in isophote twisting but that many of the galaxies stud¬ 
ied had evidence of a central spheroid. 

Irrespective of the exact origin of isophote twisting, the range 
of models able to reproduce this phenomenon is poor. The occur¬ 
rence of isophotal twists was one of the main motivations for devel¬ 
opment of triaxial Stackel models. So, it was a surprise when Franx 
(1988) showed that these models - unusually for triaxial systems 
- do not produce isophote twisting at any viewing angle. There¬ 
fore, we have hitherto been limited to constracting fully numerical 
(Schwarzschild or M2M) models to explore this phenomenon. 

Here, we demonstrate that the model explored in Section 3.1.1 
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Figure 5. Density and shape of triaxial Binney isochrone model: the left panel shows a slice of the density in the (x, y) plane and the middle panel shows 
the equivalent in the (x, z) plane. The blue colours indicate the base 10 logarithm of the density relative to the central value. The insets show zoom-ins for 
|x| < l,y < |1| and |x| < l,\z\ < 1. In the right-hand panel we show the axis ratios of the fitted ellipses in these two planes with solid lines [b/a 
corresponding to (x, y) and c/a (x, z)]. The dashed lines show the equivalent for the potential. The positions are in units of ro and the density in units of 
M/rg. 



Figure 8. Logarithmic projected density of triaxial WEH model: the model 
is viewed along a unit vector with spherical polar angles (0,6) = (j, ^). 
The lines show the major axis of ellipses fitted to the four innermost con¬ 
tours to indicate the clockwise twisting of the density contours as we move 
outwards. The positions are in units of rg and the unit of density is M/vq. 


exhibits isophote twisting. In Figure 8, we show the projected den¬ 
sity for the triaxial WEH model analysed in Section 3.1.1 when 
observed along a unit vector at spherical polar angles (0, 9) = 
{j, j). We also show the major axes of ellipses fitted to the four 
innermost density contours. Clearly, the major axis of the density 
contours twists clockwise as we move out from the centre. The gra¬ 
dient of the twist is ~ 8 ° /rg. Recalling that the effective radius of a 
Hemquist model is « 1.8rg, then the isophotal twisting gradient is 
« 14 °/ effective radius, which is very comparable to the data (see 
e.g.. Figure 8 of Leach 1981) 


4.2 Radially unstable models 

The radial orbit instability has received much attention in the litera¬ 
ture from both theoretical and numerical studies (see Merritt (1999) 
for a nice review). The underlying physics is that when spherical 
models contain a high fraction of radial orbits, the models are un¬ 
stable and tend to collapse to triaxial distributions. Polyachenko & 
Shukhman (2015) present a discussion of the mechanism by which 
the instability operates. Some systems are unstable to fast-growing 
even and odd tangential Jeans instabilities due to the lack of tan¬ 
gential kinetic energy. Other systems, which have more centrally 
concentrated radial orbits, are unstable to slow-growing bar-like 
modes of the type described by Lynden-Bell (1979) that act to align 
orbits. However, determining exactly when the instability will set 
in has been difficult to assess. Numerical experiments require the 
construction of W-body realizations so are limited to specific an¬ 
alytic equilibria such as Osipkov-Merritt f{Q) models (Merritt & 
Aguilar 1985; Dejonghe & Merritt 1988; Meza & Zamorano 1997; 
Perez et al. 1996) or polytropes / cx {—EyL'^ (Barnes et al. 
1986), or use models constructed via Schwarzschild’s method. Re¬ 
cently, Antonini et al. (2009) and Vasiliev & Athanassoula (2012) 
have showed that triaxial models constructed via the Schwarzschild 
method are also susceptible to a radial-orbit instability when the 
central anisotropy of a Dehnen model with an inner density slope 
of 1 rises to /3 « 0.4. Theoretical studies have inspected the grow¬ 
ing modes using the matrix method of Kalnajs (1977). Palmer & 
Papaloizou (1987) showed that all radially-anisotropic polytropes 
were unstable, whilst Saha (1991, 1992) and Weinberg (1991) de¬ 
rived conditions under which f{Q) models were radially unstable. 
The conclusions from both approaches do not present a uniform 
picture and it seems difficult to predict exactly when a model is 
unstable without significant computation. Our new approach, how¬ 
ever, provides another more general way to construct models that 
could be used to shed more light on the instability. 

The degree to which the radial orbit instability is important in 
the landscape of galactic formation and evolution today is unclear. 
Some galaxies are believed to be weakly triaxial but the accretion 


© 2015 RAS, MNRAS 000, 1-16 
































10 


J. L. Sanders and N. W. Evans 


Binney flattened isochrone, = 1.3, = 1.6 

0.056 



0 1 


2 3 
X 


T 


• 

poo o 

boo o 

••o o 
••• ,o 


o 

o 

o 


0.048 

0.040 

0.032 

0.024 

0.016 

0.008 

0.000 


0.232 

0.224 

0.216 

0.208 

0.200 

0.192 

0.184 


0 1 


2 3 
X 




0.056 

0.048 

0.040 

0.032 

0.024 

0.016 

0.008 

0.000 






| 0.056 

^ 0.048 

0.040 

0.032 

0.024 

0.016 

0.008 

0.000 


0.228 

0.222 

0.216 

0.210 

0.204 

0.198 

0.192 

0.186 

0.180 


0.0 0.2 0.4 0.6 0.8 1.0 

X 




0.048 

0.040 

0.032 

0.024 

0.016 

0.008 

0.000 



10 ^ 



0.0 0.2 0.4 0.6 0.8 1.0 

X 


Figure 7. Velocity dispersions and ellipses of triaxial Binney isochrone model: in the top row, the first panel shows in the [x, y) plane, the second panel 
shows Gyy in the (x, y) plane, the third panel shows in the (x, z) plane and the fourth panel shows in the {x, z) plane. In the second row, the first 
panel shows the velocity ellipses in the {x, y) plane, the second panel shows the tilt with respect to the radial direction (drawn in light grey in the first panel) in 
the (ai, y) plane, the third panel shows the velocity ellipses in the (x, z) plane and the fourth panel shows shows the tilt with respect to the radial in the {x, z) 
plane. The colours show the magnitude of the velocity dispersion along the major axis of each ellipse and only the shape of the ellipse is important. The sizes 
of the ellipses are ai'bitrary. The bottom row shows a zoom-in of the second row. The positions are in units of ro and the dispersions are in units of y/GMjr^. 


of baryons and formation of black holes may decrease the degree 
of triaxiality at the centre of galaxies. Therefore, it is conceivable 
that the radial orbit instability may be mainly of theoretical interest 
and is often damped away in realistic systems. However, in the hi¬ 
erarchical picture of galaxy formation a galaxy grows through the 
accretion of smaller satellites which can be on very eccentric ra¬ 
dial orbits. Therefore, a typical anisotropy profile is isotropic at the 
centre tending towards radial in the outskirts. It is therefore con¬ 
ceivable that the radial orbit instability sets in during the evolution 
and so is one of the mechanisms controlling the shape of galaxies. 

We limit ourselves to models where and \Jb\ are on 
an equal footing {qz = 1 in equation (16)). These should pro¬ 
duce spherical models as they only depend on Jr and L. How¬ 
ever, we allow a transition to triaxiality where the actions take 
on a slightly different meaning such that the models are not nec¬ 
essarily spherical. We construct WEH models with varying outer 
anisotropy by adjusting the parameter D \. All other parameters are 
as given in Williams & Evans (2015) isotropic Hernquist model, 
but T\ is adjusted according to equation (18). In Figure 9, we 
show the anisotropy profiles for seven models restricting the mod¬ 
els to spherical symmetry. We show the density contours and den¬ 
sity and potential shapes in Figure 10 of four of these models that 
have been given the freedom to relax to triaxiality. We see that the 


Di = 0.4684 model is spherical. D\ = 0.2 has a triaxial struc¬ 
ture in the outskirts of the density profile but the potential contours 
remain approximately spherical. For D\ = 0.114, the model has 
developed a strong prolate shape in the outskirts and the potential 
contours are prolate throughout the explored volume. When Di is 
reduced further to 0.05, the contours become more strongly pro¬ 
late. In the lower panel of Figure 9, we show the (6/a)$ shape at 
r = ro, r = lOro and r = 50ro against Di. We see that beyond 
Di = 0.2 the models become prolate in nature. The anisotropy 
profile of the spherical equivalent of the critical D\ =0.2 model is 
shown in red in Figure 9. 

It appears that when D\ > 0.2 the models remain approxi¬ 
mately spherical, whereas for Di < 0.2 the models begin relaxing 
to a prolate shape. This threshold corresponds to an anisotropy at 
the scale radius of ^(ro) ~ 0.25. We interpret this result as follows: 
for Di < 0.2 the spherical models are unstable and susceptible to 
the radial orbit instability, such that the models relax to a radially- 
distended density distribution. As we have adiabatically relaxed the 
models, we are allowing instabilities that grow adiabatically slowly 
to persist. However, we have only allowed the models to transition 
through a range of triaxial models such that, although the models 
are radially unstable, it is not clear that they will relax to the given 
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Figure 9. Radially unstable models: the top panel shows anisotropy profiles 
for seven spherical WEH models labelled by the parameter, Di, that con¬ 
trols the outer anisotropy. Spherical radius r is in units of ro. Di = 1 is 
the model presented in Williams & Evans (2015) and the other models use 
identical parameters apart from Di. The unlabelled red line corresponds to 
Di = 0.2, which is the critical model that separates the weakly spherical 
and triaxial endpoints. The bottom panel shows the shape of the potential 
contours in the {x, y) plane at three different major axis r values against 
Di for models allowed to relax to triaxiality. 


distributions. It may be that other equilibria can be explored on the 
way to the final triaxial distribution. 

In our approach the final models are approached via a series 
of triaxial non-self-consistent equilibria. However, in reality when 
D\ < 0.2 the models will pass quickly through a series of non¬ 
equilibrium states to settle down to a final, possibly quite different, 
triaxial endpoint. Therefore, we conclude that Di « 0.2 is the 
transition point between stable and unstable spherical models, and 
that for D\ ~ 0.2, we can find equilibrium models that are triaxial. 

Antonini et al. (2009) and Vasiliev & Athanassoula 
(2012) found an approximate threshold for stability of triaxial 
Dehnen models constructed with the Schwarzschild method was 
2T^/Tf ~ 1.4 where Tr is the radial kinetic energy and Tt the 
tangential kinetic energy. Assuming a constant velocity anisotropy 
this threshold corresponds to j3 ~ 0.3 (approximately the critical 
anisotropy at the scale radius) so seems consistent with the result 
found here. 

Our method involves the adiabatic relaxation of the models. 
Polyachenko & Shukhman (2015) have shown that some mod¬ 
els are susceptible to fast-growing radial modes that grow on 
timescales of order the radial period of the responsible orbits and 
Antonini et al. (2009) present models that have modes that grow on 
times of ~ 20 crossing times. As our models are generated adia- 
batically, fast-growing modes may be suppressed. 


5 DISCUSSION 

The triaxial Stackel models are exactly integrable and so possess 
no chaotic orbits. Our triaxial models are not integrable and will 
always contain some chaotic orbits, which are nonetheless assigned 
actions. In practice, the action label is meaningless and will vary 
with time. However, we have assigned the weight of this orbit based 
on the DF f{J), so the weight should correspondingly change in 
time. However, it cannot, and so the model is not in equilibrium if 
it contains chaotic orbits. 

Chaotic orbits appear regular on short times but slowly diffuse 
(Arnold diffusion). Arnold diffusion seems to be relatively slow 
compared to a Hubble time such that the equilibrium will be valid 
on scales of physical interest. However, at the centre of the models 
many of the orbits can diffuse very rapidly from one semi-regular 
orbit to the next such that the equilibrium is only valid on very 
short timescales (e.g., Valluri & Merritt 1998). Therefore, for our 
approach to hold good, we require the impact of chaotic orbits to 
be minimal otherwise the model will not be in true equilibrium. 

In fact, all equilibrium construction approaches (including 
Schwarzschild or M2M) struggle to handle chaotic orbits. In the 
case of Schwarzschild modelling, one is limited to integrating all 
orbits for a fixed period of time such that a chaotic orbit can never 
fully be incorporated. Chaotic orbits have a single integral of mo¬ 
tion, the energy, such that a model of solely the chaotic orbits can 
be constructed as fc{E). Merritt & Fridman (1996) adopt the ap¬ 
proach of summing the densities of a series of chaotic orbits in¬ 
tegrated for 100 dynamical times at each energy to construct this 
component. A full Schwarzschild model is then the combination 
of this chaotic piece plus an additional regular piece that can be 
included in the usual way. We could choose to follow the same pro¬ 
cedure here by segregating chaotic orbits and constructing f{J) 
models for the regular region of the phase-space. In this respect, 
our modelling approach does not seem inferior to other approaches 
and is probably equivalent. 

Several authors have investigated the impact of chaotic or¬ 
bits on numerically constructed Schwarzschild models. Voglis et al. 
(2002) and Capuzzo-Dolcetta et al. (2007) showed that, although 
non-rotating triaxial models contained 20 — 30 per cent of chaotic 
orbits by mass, only a few percent of these orbits by mass ex¬ 
hibited Arnold diffusion over a Hubble time. The impact of the 
chaotic orbits is a slow evolution of the axis ratios of the halos 
over time as corroborated by Zorzi & Muzzio (2012) and Vasiliev 
(2013a). However in contrast to this, Poon & Merritt (2002) have 
constructed models of triaxial nuclei with central black holes with 
considerable chaotic contributions which have fixed axis ratios 
over several dynamical times. When figure rotation is included 
the fraction of chaotic orbits by mass can increase significantly 
to ~ 50 per cent as shown by Voglis et al. (2006) and Manos & 
Athanassoula (2011). In conclusion although many orbits are for¬ 
mally chaotic in realistic non-rotating galactic potentials only a 
small fraction are chaotic on timescales of physical interest and 
these act to only weakly alter the shapes of the resulting models. 
Therefore, we conclude that the impact of chaos on the results pre¬ 
sented in this paper is weak. 


6 CONCLUSION 

The range of triaxial equilibria for stellar systems remains largely 
uninvestigated. The variety of intrinsic shapes and anisotropies for 
which stable dynamical equilibria exist is not known, and the role 
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Figure 10. Range of triaxial WEH models with variable outer anisotropy. Each row shows the results for models of differing D\ which is noted above each row. 
The left two panels in each row show density slices in the (x, y) and (x, z) planes with the insets showing zoom-ins of \x\ < 1, \y\ < 1 and |x| < l^\z\ < 1. 
The right panels show the axis ratios of ellipses fitted to the density and potential contours in these two planes [b/a corresponds to the {x, y) plane and c/a 
the (x, z) plane]. For Di > 0.2 the models are spherical, whilst for D\ < 0.2 the models collapse to triaxial distributions. The positions are in units of rp 
and the density in units of M/r^. 
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of figure rotation and tumbling is completely unexplored. One rea¬ 
son for this is the paucity of tools for building such models. 

Much of our insight into triaxial equilibria comes from the 
Stackel or separable potentials made famous by de Zeeuw (1985). 
Here, there are four orbital families, and self-consistent models can 
be built by orbit superposition or Schwarzschild’s method (Statler 
1987). There is strong evidence that the presence of four orbital 
families gives ample opportunity for exchanging orbits of differ¬ 
ent shapes, whilst keeping the total density unchanged (Hunter & 
de Zeeuw 1992). Therefore, the range of possible triaxial Stackel 
equilibria may be very rich. Other than Stackel models, there are 
only a handful of idealised triaxial equilibria in the literature (Free¬ 
man 1966; Vandervoort 1980). Of course, Schwarzschild (1979) 
and made-to measure (Syer & Tremaine 1996) techniques do in 
principle allow for a systematic exploration of triaxial equilibria. 
In practice, these techniques have been so far largely limited to ax- 
isymmetric systems because a smaller investment in computer time 
is required to sample phase-space in this case. 

Here, we have developed a new way of studying triaxial sys¬ 
tems. Many of the tools to do this exist in the literature, includ¬ 
ing algorithms for the swift calculation of actions (Binney 2012; 
Sanders & Binney 2015) and possible ansatzes for forms of the dis¬ 
tribution functions (DFs) in terms of the actions for both cored (Bin¬ 
ney 2014) and cusped (Williams & Evans 2015; Posti et al. 2015) 
galaxies. Synthesizing these results, we have shown how to build 
triaxial models with analytic DFs for the isochrone and Hernquist 
density laws. Their velocity ellipsoids are aligned with Cartesian 
coordinates near the centre, but spherical polar coordinates at large 
radii. In terms of orbital families, box orbits make a reasonable con¬ 
tribution at all radii in isochrone models, but they are less important 
at large radii in Hernquist models. 

We have studied two particular features of the models. First, 
they can exhibit isophote twisting. This occurs when the position 
angle of the major axis of the isophotes changes with radius, with 
changes of the order of 10° over an effective radius being typi¬ 
cal (Leach 1981). By projecting a triaxial Hernquist model built 
with a DF using the Williams & Evans (2015) ansatz, we demon¬ 
strated that our models could reproduce such a phenomenon. This 
is particualrly valuable, as the Stackel models - unusually for tri¬ 
axial systems - do not exhibit isophote twisting (Franx 1988). Sec¬ 
ondly, we used our models to explore the radial orbit instability. 
We considered the adiabatic relaxation of a family of DFs with cen¬ 
tral isotropy and variable outer anisotropy governed by a parameter 
Di. The endpoints of the relaxation may change from spherical to 
triaxial for some critical value of Di. Although the evolution is 
constrained, and so the endpoints may not be realised in practice, 
this nonetheless does give us the critical value at which the radial 
orbit instability sets in. We have given a particular example of this 
procedure for the triaxial Hernquist models. 

There are some further directions and possible applications 
to pursue. The study of triaxial equilibria could be advanced by 
the combination of this work with W-body simulations or pos¬ 
sibly made-to-measure models. As suggested in Binney (2014), 
these models can be used to create seeds for W-body simulations. 
These could then be evolved to investigate stability or the effects 
of growth of a disc or black hole in the model. It has been shown 
by Antonini et al. (2009) that equilibrium models constructed with 
Schwarzschild models can be dynamically unstable and it is ex¬ 
pected that some models constructed via the method presented here 
will also be unstable. The only way to address this issue is through 
A^-body modelling. Additionally, triaxiality seems an important ef¬ 
fect to account for in the interpretation of observations on external 


galaxies. To date only the Schwarzschild models of van den Bosch 
et al. (2008) have successfully modelled a triaxial elliptical galaxy. 
Therefore, one immediate application is to construct fits to ATLAS- 
3D data for slow rotators which may exhibit signatures of triaxial¬ 
ity. Finally, our models can be generalized to include streaming 
motion. The integrals of motion depend upon vf such that the sign 
of Vi is irrelevant. We are therefore able to give every particle a 
positive velocity such that the model has a net streaming velocity. 
Note that this would affect both the kinematical properties and the 
stability of the models constructed. 
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APPENDIX A: CHOICE OF COORDINATE SYSTEM 


In this appendix we detail how we choose the coordinate param¬ 
eters a and /J when using the Stackel fudge algorithm of Sanders 
& Binney (2015). Given a general triaxial potential the only inte¬ 
gral of motion we are in a position to calculate is the energy, E. 
Therefore, we choose to make the choice of a and P a function of 
E and we estimate the parameters using the closed loop orbits at 
each energy. 

We construct a logarithmically spaced grid in energy from 
Eo = cl>(0,yinin,0) to Ex = ^(Ojymax.O) with Ne = 24 grid 
points. We choose j/min = O.OlSro and i/max = ISOro where ro 
is the scale of the target model. To find the short-axis closed loop 
orbit at each energy we launch particles at a location i/i along the 
y-axis with velocity Vxi = y/'^E — 2$(0, yi, 0) in the x-direction. 
We then record the y value of the next time the orbit crosses the 
y-axis (y = yy) and the corresponding 2 ;-velocity v^f- We repeat 
this procedure and minimise 

^E(-yi - y/f + {-Vxi - Vxff, (Ai) 


where He is the frequency of a circular orbit with energy E in a 
spherical potential with the radial profile corresponding to the pro¬ 
file along the x-axis in the triaxial potential. The same procedure is 
employed to find the long-axis loop orbits by replacing x with a in 
the above equations. 

The long and short axis closed loop orhits are confined to the 
(y, z) and (x, y) planes respectively. In a Stackel potential they 
follow elliptical lines of constant A given by 


\ + P^ X + 'i 

2 2 

X y 

A J- a A J- /? 


(A2) 


With the closed orhits in a general potential found we can find the 
axis intercepts for the short axis loop: {xa,ys) and the long-axis 
loop (yi, zi). Fitting ellipses to these points we find 


P = i-y\ + zi, 
a = P- xl+yl- 


(A3) 


In Sanders & Binney (2015) we performed this procedure by 
finding the closed loop orbits around the long and short axis and 
then minimising the variation in the action describing the circula¬ 
tion of the axis as a function of the coordinate parameters. This 
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procedure was unnecessary as, with the closed orhits, found the co¬ 
ordinate parameters can be simply calculated without finding the 
actions. The procedure presented here gives identical results to that 
employed in Sanders & Binney (2015) but is less cumbersome. 

At some energies closed loop orbits cannot be found so we 
skip these points and fill them using a linear extrapolation. Addi¬ 
tionally, we constrain a < /3 < 7 such that x is the long axis 
and 2 the short axis. If our procedure produces estimates that vio¬ 
late this (e.g. in a near-spherical potential) we set a = /) — ro /10 
and /3 = 7 — rg/ 20 , where rg is the scale-length of our models. 
Decreasing this arbitrary choice of 7 — a and 7-/3 does not signif¬ 
icantly affect the accuracy of the calculations for radii r < rg but 
produces errors in the outskirts of the models (r ^ rg) as the fo¬ 
cal length is significantly smaller than the typical orbital scale. Our 
choice was selected such that the spherical isochrone and Hernquist 
models were reproduced. 


APPENDIX B: MULTIPOLE EXPANSION 

To find the potential from the density of an f{J) model we use 
a multipole expansion. Here we give details of this approach. We 
evaluate the density on a 3D grid in spherical polar coordinates 
(r, (j), cos 9). The grid in radius r is logarithmically spaced between 
two radii rmin and rmax such that 

ri = aoe^' -f ^ min 5 (Bl) 

where 

Si = jpZTi ^max/^min); "— 0, 1 . (B2) 

The scale ag is chosen to be the scale radius of the model of interest, 
rg, and we set rmin = O.Olrg and rmax = 200rg which encloses 
99 per cent of the total mass in both the spherical Hernquist and 
isochrone models. The cos 9 and (j) grids are based on Wa-point 
Gauss-Legendre quadrature on the intervals (—1,0) and ( 0 , 7 t/ 2 ) 
respectively. We set Nr = 40 and Na = 15. 

The potential is expanded as 

^ —^max 771=1 

<E>(r,<(.,0) = - 47 tG ^ ^ 0,m(r)yr(<^,0). (B3) 

!=0,2,... m=0,2,... 

where Yi^{<j),9) are spherical harmonics. The sum is over even 
values of I and m using a maximum I of Imax = 8. We choose 
to work with a symmetric definition of 9) (Vasiliev 2013b) 

given by 

Yr'{(t>,9) = ^ Pr (cos 9) cos mcf), (B4) 

where y = v/2 if m > 0 otherwise (V = 1. The normalized 
associated Legendre polynomials, P™ (cos 9), are calculated using 
the GNU Science Library (GSL, Galassi et al. 2009) as 

l £i\ / 2Z A 1~ / (/ ?7l)! 

where P™(cos0) are the associated Legendre polynomials. The 
Y{"{(f), 9) satisfy the orthogonality condition 

f d9 f d(t) sm 9Yr'Y iY' = (B 6 ) 

Jo Jo 2 t +1 

4>im(r) is found by integrating the potential of the series of 


spherical shells of thickness Sr with mass pim{r)Sr inside and out¬ 
side the radius r. Mathematically we have 

<(>i™(r) = r“'“^7™(r) + (B7) 

where 


daa+^pim(a), 

^0 

/ CX5 

daa"'+Vim(a), 


(B 8 ) 


Pim(r) = 8(-l)""/^y d(cos9) j d(f)Yr{(j),9)p{r,(j),9). 


(B9) 


The integral for (pim (r) is calculated using the trapezoidal rule in 
the coordinate S = log((ri — rmin)/ao) and the double integral for 
pim{r) is found using Gauss-Legendre quadrature. 

Similarly, to find the forces we must evaluate 


dcplm 

dr 


1 + 1 


f) + Ir 


1 — 1 t(oo) 
^Im 


f), 


(BIO) 


and the derivatives with respect to (f) and 9 are found by analytically 
differentiating Y/^{(f>, 9). 

We begin the calculation of these quantities by evaluating the 
density on the grid in (r, (p, 9). This is the costliest part of the calcu¬ 
lation but it can be parallelized. Then for each I and m we evaluate 
pim{r) on the grid in r and use these results to compute (f>im{r) and 
dpim/dr on a grid in r, I and m. Assuming a finite mass model (i.e. 
pim{r) has a shallower divergence than as r —>■ 0 and falls off 
steeper than r~^ as r —>■ 00 ) the integral from r = 0 to r = rmin 
is found as I,™ (rmin) = |pim(j'min)fmin- We Set 7 /“^ = 0 for 
r A rmax. 

For r < rmin we use quadratic extrapolation of (j>im{r) and 
linear extrapolation of dpim/dr. For r > rmax we assume the 
density is zero and so extrapolate (j>im{r) and dpim/dr as 


+ t'max) — 0im(fmax 

) = 6, (r 
or 



(BID 


The potential and forces calculation was tested using a triaxial 
perfect ellipsoid Stackel potential (de Zeeuw 1985) as well as a flat¬ 
tened Hernquist model (replacing the square of the spherical radius 
r^ with m? = + {y/qyf + {zjqzf in the Hernquist density 

profile) for which the potential and forces were calculated using 
equation (2.140) of Binney & Tremaine (2008) for the potential of 
a general density distribution stratified on concentric ellipsoids. 


APPENDIX C: CROSS-CHECK WITH GENERATING 
EUNCTION APPROACH 

In the main section of the paper we have used the triaxial Stackel 
fudge approach of Sanders & Binney (2015) to construct the mod¬ 
els due to its speed. However, the speed comes at the cost of in¬ 
creased errors in the actions. Using the Stackel fudge the action 
estimates will oscillate. For some orbits this oscillation is genuine 
(see Section 5 on resonant and chaotic orbits), whilst for others it 
is a natural consequence of the approximation used. This naturally 
will introduce errors in the calculations to construct the models. 
Sanders & Binney (2015) explored this and discussed the error the 
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Stackel fudge method produced when finding the moments of simi¬ 
lar non-self-consistent models. It was found that the Jeans equation 
was satisfied to < 10 per cent for two tracer DFs in a triaxial NFW 
potential. 

Sanders & Binney (2014) presented an alternate method for 
estimating actions in triaxial potentials. This algorithm will con¬ 
verge on the true action when one exists (for regular orbits), but is 
unfortunately much slower as it relies on orbit integration. However 
it can be used here as a cross-checking tool. Here we summarise the 
method and present a calculation of the density of one of our mod¬ 
els using this approach. 

Cl Generating function construction 

The generating function method of Sanders & Binney (2014) is 
based on the construction of orbital tori methods pioneered by 
McGill & Binney (1990) and operates by construction of a gener¬ 
ating function from the angle-actions in a simple analytic potential 
to those in the target potential from a series of orbit samples. Here 
we very briefly detail the method. 

Given a target potential and a initial {x,v) coordinate we in¬ 
tegrate the orbit. Here we choose to integrate for T = Sforb where 
forb is the orbital period of a circular orbit with the same orbital en¬ 
ergy in a spherical potential with a radial profile corresponding to 
the a;-axis profile of the triaxial potential. From this orbit we sample 
Nt = 200 points and infer the orbital type (loop or box) from the 
angular momentum. We then choose a corresponding toy potential 
(triaxial harmonic oscillator for the box orbits and isochrone for 
the loops). We estimate the parameters of the toy potentials by fit¬ 
ting the toy forces to the true forces at the maximum (|x|,|t/|,| 2 |) 
points for the harmonic oscillator and by fitting the toy forces to 
the true forces at pericentre and apocentre for the isochrone poten¬ 
tial (note this differs from the procedure employed in Sanders & 
Binney (2014) who minimised the variance of the toy Hamiltonian 
for the orbit samples). With the toy potential estimated we find the 
toy actions and angles {J' ,G'). The true actions are then related to 
these toy variables via the generating function S{9', J) as 

J — j' — 2 nSn cos n ■ 9', (Cl) 

n 

where Sn are unknown Fourier components of the generating func¬ 
tion. From our orbit samples we have Nt such equations which can 
be solved for {J, Sn) by minimising the sum of the square differ¬ 
ences between the toy actions and those calculated using the gen¬ 
erating function series. We consider terms up to |n| ^ A^max = 6. 
In addition to the checks of the solutions employed in Sanders & 
Binney (2014) to ensure good coverage of all modes, we also en¬ 
sure the radial action for the loop orbits is positive, all actions for 
the box orbits are positive and the average of the action over the 
corresponding angle does not deviate from the estimated action by 
more than the maximum estimated action. If these criteria are not 
satisfied we increase the number of samples (if there is less than 
one sample per tt phase of n • 0' for any mode n) or the time inte¬ 
gration window (if n • 9' does not wrap at least 27 r for any mode n) 
by a factor of two until an upper limit (four times the initial choice) 
and if the conditions are still not satisfied we return the average 
actions. 



10'^ 10-' 10“ lo' 
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Figure Cl. Density profile of triaxial WEH model along a ray with spheri¬ 
cal polar angles (</>, 6) = (^, ^) calculated using the triaxial Stackel fudge 
method (black solid) and the generating function method (red dashed). The 
lower panel shows the relative difference in the density. The positions are 
in units of tq and the unit of density is M/r^. 

Figure C1 we plot the density profile of the converged WEH model 
of Section 3.1.1 along a ray at spherical polar angle (4> = j,8 = 
calculated using the fudge method and the generating function 
method, and in Figure C2 we show the same calculation for the 
isochrone model of Section 3.2.1. Note the full convergence pro¬ 
cess was not performed using both approaches. We simply took 
the converged potential from the Stackel fudge calculation and re¬ 
calculated the density in this potential using the generating func¬ 
tion approach. For the WEH model the relative error in the den¬ 
sity is < 10 per cent everywhere. The largest errors occur outside 
r « 5ro. For the isochrone model the relative error in the density is 
< 2 per cent for r < 5ro but the error steeply rises beyond this up 
to 30 per cent for r = 20ro. It is unclear whether this systematic 
discrepancy is due to the fudge or generating function calculation 
but it appears to occur in the regions where the long-axis loop or¬ 
bits are prominent so perhaps is due to their mishandling. We note 
that inspecting a less extremal model (a^f, = 1.3, Oz = 1.4) we 
find that the two methods are agree well. 


C2 A cross-check 

Here we perform a cross-checking calculation for two of our mod¬ 
els to measure the error the Stackel fudge approach introduces. In 
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Figure C2. Density profile of triaxial isochrone model along a ray with 
spherical polar angles (0,0) = (J,^) calculated using the triaxial 
Stackel fudge method (black solid) and the generating function method (red 
dashed). The lower panel shows the relative difference in the density. The 
positions are in units of rp and the unit of density is Mjr^. 
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